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We present a new dynamical model of DNA. This model has two degrees of freedom per base- 
pair: one radial variable related to the opening of the hydrogen bonds and an angular one related 
... to the twisting of each base-pair responsible for the helicoidal structure of the molecule. The small 

' amplitude dynamics of the model is studied analytically : we derive small amplitude envelope 

solutions made of a breather in the radial variables combined with a kink in the angular variables, 
showing the role of the topological constraints associated to the helicoidal geometry. We check the 
stability of the solutions by numerical integration of the motion equations. 
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I. INTRODUCTION. 



Simple dynamical DNA models have been recently studied ^ to understand the possible basic mechanisms of some 
biological processes such as denaturation, initiation of the transcription, and transcription. The thermal denaturation 
has been investigated by Peyrard, Bishop and co-workers with a unidimensional model (PB model) which allows 
' local openings of the hydrogen bonds and formation of denaturation bubbles. The local openings can be analytically 
described as breather-like objects of small amplitude, which have nevertheless interesting properties: as long as their 
. amplitude is small enough, they can move along the chain, collect energy and grow They can also be trapped 
' by some local dishomogeneities [D , which suggests that the properties of breathers could allow the formation of the 
. transcription bubble after the interaction with the bound RNA-polymerase. 

However, in a realistic description of DNA dynamics, the topological constraints related to the helicoidal structure of 
the molecule cannot be ignored. Activation or repression effects caused by conformational changes and, in particular, 
^ ■ prevention of transcription due to a large positive excess in twist, are known and largely investigated in biology. In 
fact, during all processes in which the DNA base-pairs open, a local unwinding of the helix follows for topological 
reasons. Consequently a local extra-twist accumulates at the two ends of the bubble and induces a long range elastic 
Q ' stress. Although the mechanical properties of DNA have recently been the object of a renewed interest (6|j^, previous 
, attempts to take into account the helicoidal structure in the models have been limited to the introduction of the forces 
' that can appear due to the proximity in space of bases which are non-adjacent in the sequence |^]. 
\ Our aim here is to build a simple model that takes into account the twist-opening interactions c ue to the helicoidal 
^ ' molecular geometry. Such a model provides an extension of the PB approach towards a more realistic description of 
k> biological processes. It can also be useful for more general studies of the interaction between geometrical conformation 
, and dynamical properties of the molecule referring to the recent mechanical experiments on DNA |8|-[l0[| . 
' In designing such a model each base will be considered as a single, non-deformable object. We introduce two 
degrees of freedom per base-pair: one radial variable related to the motion of the bases along the diameter that joins 
their attachment points to the helicoidal backbone, and a twist angle of each base-pair defined by the angle between 
this diameter and a reference direction. This twist angle, which increases from one base-pair to the next one, is 
responsible for the helicoidal structure of the molecule. From the resulting Lagrangian the equation of motion for 
small displacements with respect to the equilibrium position are written. We then derive analytical approximated 
solutions for the small amplitude nonlinear distortions of the molecule. As suggested by the helicoidal geometry we 
show that the radial breather-like opening is associated to an angular un-twisting of the molecule. We finally present 
the results of numerical simulations of the dynamics showing that these solutions are stable for long periods of time 
and can move along the chain. 

The analytical method of derivation of the envelope soliton solution of the equation of motion used in this work 
has been proposed and developed extensively in a previous paper ]l5| . The reader is referred to this work to better 
follow the technical steps of our calculation. 
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II. THE HELICOIDAL MODEL. 



In the PB model , the bases are point masses allowed to move only in the direction of the hydrogen bonds that 
connect them. A Morse potential describes the effect of these bonds, while neighboring bases along the same strand 
are harmonically coupled, simulating the stacking interactions. 

Similarly, in our model the group made of a sugar ring and its connected base is treated simply as a point mass 
(without distinction between the different base types); the phosphate backbone between two base-pairs is modelised 
as an elastic rod. The additional twist motion is now introduced by allowing the two bases in each pair to move in 
the base-pair plane instead of constraining them on a line. It is convenient to choose a polar coordinate system (Fig. 
1). The model does not attempt to describe the acoustic motions of the molecule since only the stretching of the 
base-pair distance is considered. This amounts to fixing the center of mass of the base pair, i.e. the two bases in a pair 
are constrained to move symmetrically with respect to the axis of the molecule. Then, to describe the stretching of a 
base pair and the variation of the helicoidal twist we need only two degrees of freedom per base pair: the coordinates 
r„ and Lpn of one of the two bases with respect to a fixed reference frame. 

As in the PB model, a Morse potential describes the hydrogen bonds linking bases in a pair with an equilibrium 
distance i?o- 

A proper choice of the coupling between radial and angular variables has to reproduce the equilibrium helicoidal 
structure. In DNA, the latter originates from the competition between the hydrophobic effect (that tends to eliminate 
water from the core of the molecule by bringing the neighboring base-pair planes closer) the electrostatic repulsion 
between neighboring base planes (which has the opposite effect) and the rigidity of the two strands (that separates 
the external ends of the base pairs by essentially a fixed length related to the phosphate length) jl^ . We combine the 
first two forces in a unique stacking effect that fixes base-plane distance h in the model. As the equilibrium backbone 
length L is greater than ft,, it is then necessary to incline the strands in the typical helicoidal structure to minimize 
the energy. 

The geometrical parameters h and i?o have been chosen to match the structure of B-DNA. Then, in order to impose 
an equilibrium twist angle — ^n-i — i0o between two consecutive base pairs, we select the equilibrium length L 
of the springs representing the phosphate backbone to be 



L ^ ^ + iRlain^-^) > h (1) 



However, because of the invariance of L with respect to the sign of Oq, the elastic rods rigidity is not sufficient to 
guarantee the correct helicoidal shape. A zig-zag structure with a random succession of ±0o would as well minimize 
the energy of the system. This can be avoided by adding a three-body curvature term to the Lagrangian, that imposes 
a continuity between the differences ifn — ^n-i and ipn+i — ^n- 

The final Lagrangian is written considering the expression of the three-dimensional distance between the two bases 
along the strand (Fig. 1); it reads 



- Y ft2 _^ ^2_^ ^ ^2 _ 2r„_ir„ cos(i^„ - ^Pn-l) - L) 

n 

- ^ Go (</3„+l + ^n-l - 2<p„)^ (2) 

n 

where m is the base mass, D and a are the depth and width of the Morse potential well (without distinction between 
double and triple hydrogen bonds), K is the backbone elastic constant and Go the backbone curvature constant. 

The results presented below have been obtained with, m = 300m. m. a., D = QMeV, a = 4.45 A^^ that are the 
same valued adopted in the PB model, and with K — l.QeV and Go = KRq/2. We intend to refine these 
values in future statistical mechanics studies by comparing the predictions of the present model with some available 
experimental data, eg. the temperature of DNA denaturation and the rigidity of the molecule fl^ . 

For the geometrical parameters we adopt the B-DNA values Rq ~ 10.0 A, 8o — 36° and h — 3.4 A. 

One can simplify the Lagrangian by introducing the adimensional variables r'^^ — ar^, the rescaled time t' — 

Da^ /m t and the renormalized parameters 
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K' = K/Da"^ 

g' - Go/D 
R'q = chRq 

h' = ah 

L' = aL. (3) 
In the following the new variables will be written without primes. 

We perform an expansion of the energy of the backbone springs in up to the second order around the equilibrium 
position: 

Un = Tn - Ro (4) 
(l)n ^ RoiVn - nQo) , (5) 



obtaining 



/i^ + rl_-^ +rl- 2r„_ir„ cos((p„ - (fn-i) ~ L) 



R^ 2 

-|[sineo(0„ - «!>n-i) + (yn + yn-i)(l - coseo)] . (6) 



We also expand the Morse potential up to the fourth order in j/„ 

(e-^" - 1)' = \vl \vl + Y^vt + 0(y„^) . (7) 

The difference in the order of the two expansions is consistent with our parameter choice. 
With these expansions, the equations of motion become: 

-:^(0„+i-0„„i) (8) 

<Pn = i^00(0n+l + 0ri-l " 2(/)„) 

- G{4>n+2 + 4>n-2 " Hn+l - 40„-l + 60„) 

K 2 2 2 1 

H ^^(^"+1 - Vn-l) - -^VnCf^n " -^VuK - -^VuVnCt^n - -^vIK (9) 

z rto -fCo j£q rC, 



where 



i^yy = (KRl/L^) (1 - cos Oo)' , (10) 

K^^ ^ {KRl/L^) {sin^eo) , (11) 

i^j,^ = 2(/fi?2/L2)(sineo)(l -cosOo) (12) 

G - g/Rl . (13) 

The constants Kyy,K^^ are the effective elastic constants respectively for the base pairs opening and the twist 
rotation, Ky^j, is the coupling constant between stretching and twist; these three constant are geometrically related. 

III. LOCALIZED BREATHER-LIKE SOLUTIONS. 

We now look for nonlinear localized solutions, characterized by the propagation of a coherent collective structure 
on a time scale greater than the time scale of the vibrations of each particle around its equilibrium position. 

According to the method developped in we first solve for a wave packet solution of the linearized system with 
weak dispersion. This amounts to solving 
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{J{q)-Ujf{q))Vi{q) = 



(14) 



where I — +, — is the branch index and 



1 + 2Kyy ( 1 + cos q) iKyfjj sin q 



\ -iKy^sinq 2K^^{1 ~ cosq) + G{6 - Scosq + 2cos2q) J - \c* b' ^ ' 

finding the eigenvalues and eigenvectors 



col{q)^-ia + b±^ia-br+A\c\^) (16) 



V^{q)^J^±ya-^ j (17) 

where Af± are the vector norms. 

Then we apply a perturbative expansion of the system ( [1^ ) around one normal mode. We have chosen as carrier 
wave a mode {qQ,LU^ = w_|_(go)) on the optical branch and with a small wave number that corresponds to a prevalently 
radial excitation weakly oscillating. In fact, keeping in mind DNA opening, we shall calculate the behaviour of a 
localized bubble-like radial distorsion and the induced angular dynamics. We obtain the wave-packet group velocity 
from the first order eigenvalue correction 

^^+^ = ^^Z'J^ = 2^(^'\^i^\' + + c'*V+*V+ + b'\V+\') , (18) 

where the primes indicate the derivatives with respect to q calculated at qo. We then obtain the correction to the 
eigenvector 

(19) 

V-*J'V+ a'VrVi^ + c'Vi~*V2^ + c'*V,~*V+ + b%-*V+ 
a = — — = 2 2 (.20) 

W+ — 071 ^~ 

and from the second order eigenvalue correction the group velocity dispersion is obtained as 
- ~2^ + u^l-ul > 

= —ilia"\V+\' + c"V+*V+ + c"*V+*V+ + b"\V+\') + |ap(c.2 - )) (21) 

We now take into account the nonlinearity. We look for a small amplitude solution. The expansion parameter e 
for the solution is introduced to solve the equation of motion at increasing order of accuracy inserting the nonlinear 
terms in a progressive way. This iterative expansion is combined with the expansion in multiple scales for the weak 
dispersive wave packet. The latter requires the introduction of the variables xi = ex,ti = ei,t2 = for the slowly 
varying amplitudes, where e is the same parameter as before. We look then for a solution of the form 

+ e2e2'(««"°-'^+*«)7(xi, ii, i2) + e'fiixiMM) (22) 

where the first term is the wave packet, the second term arise because the linear system admits a constant solution 
with a non zero second component, corresponding to the null column of Jifi). The constant and second harmonic 
terms are induced by the quadratic nonlinearity in the equation of motion. 

We have now to determine the slowly variyng amplitudes yl(xi, ti, ^2), f?(a;i, ii, i2)j 7(2^1, ^i, ^2)1 m(2;i, ti, ^2) 
From the O(e^) equation of motion we obtain for the second harmonic terms the system of equations: 

(J(2,o) - M)7 = i '^'li' \ A- . (23) 
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Its solution may be written as 7 = jcA"^ and for the constant terms 

J(0)M^-J'(0)A,^=(3lVP+^^|V^/P) 1^1^ . (24) 

Inserting the value of J(0) from ( ]l5| ) in the above matricial equation, we obtain the identity 

(1 + AKyy)^^, + Ky^-^^a^ = (3|F+l2 + ^|V^2+n 1^1' (25) 

that does not allow to determine the two unknown variables. We have thus to consider the O(e^) system of equations 
for the constant terms 



Pi f)2 f "^1^+ \Y+\2( A* dA OA' A\ 

ij>(0)JLn- j"(0)—a - ' ^ ' ^ ' 1 (26) 



corresponding to the two equations: 



d 2iuj+ dA OA* 



The equation (|2^) gives /X2 as a function of A. The equation (^8|) can be integrated using the wave packet property 
= — cj^^^-* We then obtain from (|25|), ( p8|) the following system for /ii and 172, 



1 + 4A„ AW \ ^p, \ / 3|V+|» + ^|V+I» 



the solutions of which may be written as 

m=Micl^P (30) 
(j2^(Tc j \A\'^dxi (31) 

The nonlinear O(e^) terms in exp i{qono — uj+to) in the equation of motion give rise to the nonlinearity that balances 
the group velocity dispersion of the wave packet. We indeed obtain the following Non Linear Schrodinger (NLS) 
equation for the envelope A expressed in a frame moving at velocity to^^^ (variables S = xi ~ w^^^'ti, t — 12)'- 



with 



iAr + PAss + Q\A\^A = Q (32) 



^(2) 

Q=mr±Q^^ (34) 

2(jj+ 



where 



2 Ha 



4^2 

3^i+*7ic + -^^^2+*72, + 3V+^i,, (35) 

h = 1^(4^ - 2^2+ 71c) + ^^2^Ailc - 

^VVr + ^(2|l^r+|V+ + l^+Vr) . (36) 
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If PQ > 0, equation (|32|) has an envelope soliton solution 



1 u 
A{S, t) = A sech[— (5 — Uct)] exp [iTr^iS — u^t) 



where 



2PQ 
2P 



(37) 

(38) 
(39) 



are respectively the amplitude and the width of the curve. 

Once the NLS equation is solved for ^1(5*, r), we have a2, 7, ^ (solving ( p7| ) we obtain ^2 — M2c/ l^Pc^a;!) and then 
the complete solution 



y = e{V+ - ieVl^^^)Ae'' + c.c. + e^ji.A^e^'" + c.c. + e^id^l" + ^(e^) 

.(1) 9 



(t> = ea^ j \A\^dxi + e{V+ - ieV^'^^ -^)Ae'^ + c.c. + e^^Sc^^e^** + c.c. + eV2c ^ l^l^da^i + 0(6^) . 
where — (qona — i^+to) ■ The final result can be rewritten in the following form: 



(40) 
(41) 



y = 2eVi+ Asech[r]{x ~ Vet)] cos {ICx - fit) 

'2 



+ e2Vi(i)^sech[r;(a; - VJ)] 



tanh {r]{x — Vet)) sin (ICx — fit) + — cos {ICx — fit) 
L„ P 



2e'^jieA'^sech^lT]{x - Vet)] cos (2/Cx - 2^) 
fiiee'^A^sech^[r]{x - Vet)] + O(e^) 



(42) 



e{<Te + e^i2c)LeA^ tanh {ri{x ~ Vet)) ~ 2e^|V2+|sech[7/(a:: - Vet)] sin (/Cx - fit) 
eV2^^^y^sech[77(a; - Vet)] ^ tanh (77(0; - Vet)) cos (/Ca; - m) + ^ sin {Kx - 
2e2|72j^^sech2[ry(a; - Kt)] [sin {2ICx - 2flt)] + O(e^) 



(43) 
where 

(44) 

(45) 
(46) 
(47) 

(48) 

Fig. 3 shows the analytical solution (|42|),([43|) in the original variables (not renormalized), as a function of x and t. 
Fig. 4 displays the results of numerical simulations of the system described by the complete Lagrangian (^) starting 
from the analytical solution (^2|),(^) for f = as an initial condition. We have chosen the carrier wavevector go = 0.1. 

The shape of the solution is best seen on Fig. 3a, 3b which show for a few periods the radial displacement r„ (t) — Rq 
(Fig. 3a) and of the angular displacement (pn{t) (Fig. 3b), while Fig. 3c represents the twist angle, i.e. the difference 
between neighboring angles, A(pn{t). As in the PB model the radial motion has the shape of an asymmetric breather, 
the stretching of the base-pair distance being larger than its compression, as expected from the asymmetry of the 
Morse potential. The corresponding angular motion clearly exhibits a kink structure due to the non oscillating term 
a{x,t); a small oscillating part due to the other terms in (B3) is superimposed on it. 

The shape of the solution agrees with the geometrical properties of the molecule related to its helicoidal structure. 
It shows, as expected, that our model can indeed describe the local untwisting which should be coupled with the radial 
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breather due to the geometrical constraints. To make this clearer we plot in Fig. 3c the evolution of the twist angle. 
In addition to the local opening, small periodic overtwists on the boundaries, coming from the oscillating part of the 
solutions are visible. The coupled breather and kink solutions move with a peak velocity Ve — lo'^l' + eug without 
changing their internal structure. 

The long time evolution of the initial condition, calculated with the full Lagrangian and shown in Fig. 4, attests 
that the local excitation is very stable. It does not radiates energy with the exception of a small initial transient that 
could be eliminated with absorbing boundary conditions. During all the simulation time the total energy lost with 
the absorbing conditions is less than 10^'°^ of the total energy. It should be noticed that Fig. 4 does not show the 
actual oscillations due to the large sampling time used to produce the graph. The slow oscillation visible on Fig. 4 is 
due to a beating between the actual oscillation and the sampling time. Fig. 4 shows in fact about 2230 oscillations 
which proves the exceptional stability of the breather in the helicoidal structure and the ability of the geometrical 
nonlinearity and of the nonlinearity coming from the Morse potential to keep the energy localized. 

As our solution has been obtained with a small amplitude expansion, we should not expect it to be a good solution 
when its amplitude increases to much. This can be seen on the results of a numerical simulation starting from a much 
larger amplitude (Fig. 5). In this case the initial condition decays very fast and radiates energy before stabilizing 
to a breather which has still a pretty high amplitude and which is localized on about 6 sites (width at half height). 
Although it is very narrow the breather is still mobile in the system and it is accompanied by a sharp kink as expected 
from the geometrical constraints. The solution reached after the transient is stable as shown in Fig. 5c: although 
the simulation has been run with absorbing boundary conditions, the energy tends to stabilize once the radiated 
waves have been dissipated at the ends of the molecule. In all the solutions described above the breather frequency 
ri remains in the gap between the two branches so that the breather does not emit acoustic phonons. 

IV. CONCLUSION 

In conclusion we have proposed a two-variable helicoidal DNA model that, from a biological point of view, can 
be a better starting point than the PB model for the studies of all the DNA dynamical processes involving internal 
openings. Our model also describes the untwisting that accompanies the base-pair stretching because of the helicoidal 
molecular geometry. We analytically derive, using the method introduced in p5| , a breather like solution in the radial 
variable coupled with an angular solution dominated by a non oscillating kink-like term. A further step will be the 
study of the formation of this kind of excitation in presence of a thermal bath and of the onset of thermal denaturation. 
This will be useful in determining the model parameters by comparison with denaturation temperature experimental 
data. Furthermore we expect to find interesting thermodynamic properties due to long range effect generated by the 
helicoidal DNA structure. 

As far as biological processes are concerned, the model may however still lack some relevant features, in particular 
the possibility of bending and stretching of the molecule. In order to obtain a stretchable molecule axis we can simply 
replace the rigid constraint of fixed base planes distance h with an elastic term which leads to the introduction of a 
third degree of freedom with a new elastic constant. This could be a link to some recent experimental studies of the 
mechanical response of a DNA molecule to some twisting and stretching forces . Anyway the model contains in 

its possible dynamical behaviors some other important features, as e.g. the possibility of local transitions from right 
hand B-DNA to left hand Z-DNA, related to the curvature term in the Lagrangian depending on the constant Gq. 
Our model is, finally, a first attempt to insert in simple DNA dynamical model topological constraints that are of 
great interest in today's biological studies. 
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FIG. 1. Fig. 1: Choice of variables (left) and the schematic view (right) of the model. The n-th base-pair is represented 
with its radial and angular displacements (r„,(p„) with respect to a fixed external reference frame. Neighbouring base-pairs 
are linked along the two strands by elastic rods with equilibrium length L. Ro, Go and h are the geometrical parameters of the 
model, with values deduced from the real B-DNA geometry. 
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FIG. 3. Fig. 3a: Analytical solution (|42|), ([l3|): we show here the radial displacement r„{t) — Rg, obtained with the model 
parameters indicated below and with solution parameters e — 0.1, ue = —0.1, uc = 0.1. The typical breather oscillation is 
clearly visible; its period is of 105.1 time units {It.u. = 1.0210~^'*s). The amplitude 2eA/a and the half heigth width, which is 
about 2.77Le/e, are respectively 0.016 Aand 20b.p. 
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FIG. 4. Fig. 3b: Analytical solution (M), ([l3|): the corresponding angular displacement ipn{t). The main contribution to 
this solution is the kink-like term eacLc A tanh (rjlx — Vet)) of 
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FIG. 5. Fig. 3c: Analytical solution (^^, (^3|): the twist angle, given by the difference between neighboring angles, Aip„{t). 
We stress that the breather opening mode in the radial variables corresponds to an untwisted region. 
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FIG. 6. Fig. 4a; Numerical integration (obtained with a Runge Kutta algorithm) of the initial condition which corresponds 
to the analytical solution of Fig. 3. The amplitude and width parameters preserve their initial values 2eA ~ 0.016 Aand 
2.77Le/e ~ 20 6.p. with good accuracy, as well as the oscillation period. The total simulation time corresponds to 2234 
oscillation periods. Wc introduce hero absorbing boundary conditions to avoid noise effects produced by reflection of the small 
amount of radiated energy (less then 10~^ of the total energy). 
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FIG. 7. Fig. 4b: Angular displacement ipn{t) generated by the numerical integration of the same initial conditions. 
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FIG. 8. Fig. 4c: Numerical integration of the same initial conditions, twist 
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FIG. 9. Fig. 5a: Numerical integration of an initial condition obtained with solution parameters e — 0.2, ue — —0.2, uc — 0.6. 
We perform again the simulation with absorbing boundary conditions. The initial amplitude is 2eA = 0.1 Aand the half height 
width 2.77Le/e ~ 4 6.p. This initial solution is not stable. The initial radiation is clearly visible. However, after some time, the 
system tends to stabilize to a smaller brcatlior+kink mode with amplitude 2cA ~ 0.07 and width 2.77 Lc/e ^ 6b.p. (Wc stress 
that the apparent slow modulation of the maximum amplitudes and the apparent oscillation period are due to the beating 
between the effective oscillation and the sampling time.) 
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FIG. 10. Fig. 5b: Numerical integration for the same parameters as in Fig. 5a, angular displacement (finit). 
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FIG. 11. Fig. 5c: Total system energy during the simulation with parameters as in Fig. 5a. The loss of energy is due to the 
absorbtion at the boundaries and represent essentially the radiative part. After the initial transition to the smaller amplitude 
breather +kink mode the energy tends to stabilize. 
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